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Abstract 

Many physical, chemical and biological systems exhibit a cooperative 
or sigmoidal response with respect to the input. In biochemistry, such 
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behavior is called an allosteric effect. Here we demonstrate that a 


system with such properties can be used to discriminate the amplitude 
or frequency of an external periodic perturbation or input. Numerical 
simulations performed for a model sigmoidal kinetics illustrate that 
there exists a narrow range of frequencies and amplitudes within which 
the system evolves toward significantly different states. Therefore, 
observation of system evolution should provide information about the 
characteristics of the perturbation. The discrimination properties for 
periodic perturbation are generic. They can be observed in various 
dynamical systems and for different types of periodic perturbation. 


1 Introduction 

Bistability and hysteresis are commonly observed in physics, chemistry and 
biology IDEl. Let us assume that a system has two stable states Si and S2 
and that an increase in the value of control parameter A above the threshold 
Ai triggers the transition from Si to S'2, whereas the reverse transition from 
S2 to Si occurs if the value of the control parameter drops below A2. Such 
a system can obviously be used as a discriminator of the control parameter 
value. For example, if the initial state is Si and after some time we observe 
the system in S2, then at some point the value of the control parameter 
necessarily exceeded Ai. However, if only time-monotonic changes in the 
value of the control parameter are considered, then the system discrimination 
ability is reduced to just two values ; Ai and A2. 

In this paper we demonstrate the suitability of a dynamical system charac¬ 
terized by sigmoidal kinetics for discrimination-oriented applications, under 
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a new strategy of imposing a periodic perturbation or input on a cooper¬ 
ative system. It has been reported P-[S] that periodic perturbations can 
significantly change the time evolution of a nonlinear system. As a discrim¬ 
inator prototype, we consider a two-variable system in which the inflow of 
one of the variables is a control parameter. Numerical simulations reveal a 
non-trivial property of such a system: a marginal change in the inflow pa¬ 
rameters (amplitude or frequency) can switch the response of a cooperative 
system between different branches in the stage diagram. The frequency at 
which such switching occurs is a monotonic function of the inflow ampli¬ 
tude. Therefore, at a hxed amplitude of periodic inflow, the observation of a 
transition between different types of oscillatory evolution of the system pro¬ 
vides information which allows us to discriminate the inflow frequency. The 
above discussion does not necessarily limit the range of frequencies that can 
be discriminated by the observation of transitions between different types 
of oscillations Similarly, for a fixed frequency of periodic inflow, transitions 
between different types of system oscillations occur within a narrow range of 
amplitudes. The transition can be used to discriminate the inflow amplitude, 
but for the model considered here, the useful range of such discrimination is 
rather limited. 

In numerical simulations, we consider simple system dynamics dehned by 
a single sigmoidal term expressed by a rational function, which is typical for 
enzymatic reactions [ 9 HI 3 ]. In such reactions the appearance of sigmoidal 
kinetic behavior is usually interpreted to be the result of the interaction of 
substrates with enzymes through positive cooperative binding. Modeling of 
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cooperative binding leads to the Hill eqnation [5]: 

Ki + |L]” 

where 9 is the fraction of ligand binding sites hlled, [L] is the ligand con¬ 
centration, Kd is the apparent dissociation constant derived from the mass 
action law, and n is the Hill coefficient which represents the degree of coop- 
erativity. If n = 1, there is no cooperativity; for n > 1, the cooperativity 
is positive. Kinetics with sigmoidal behavior are not limited to enzymatic 
reactions. This also describes the response of various biological systems to 
external stimuli, including the effect of drug delivery, which is an interesting 
topic in pharmacology. Among the many experimental studies that have re¬ 
ported sigmoidal behavior, the Hill coefficient n usually has a value between 2 
and 4 [HHIS]. Here we selected n = 3 for the numerical simulations presented 
below. 

The paper is organized as follows. In the next section, we consider a 
bistable model and study its time evolution as a function of the amplitude and 
frequency of periodic inflow. We demonstrate how the system can be used 
as a discriminator and discuss the sources of discrimination errors. In the 
hnal section we argue that the observed phenomenon is generic and discuss 
its potential applications. 
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2 The response of a model dynamical system 


to periodic perturbations 


Let us consider a dynamical system of two variables (x(t), y{t)) defined by a 
set of differential equations: 
dx 


dt 


= g{x, y, t) = -ax + y + A- (sin(27r/t + 4>o) + 1) • © W (2) 


$ = h{x,y) = - ■ 


X 


-y) 


(3) 


dt e + 

In Eq.(|2]), the last term I(t) = A ■ (sin(27r/t + 0o) + 1) • 0(^) describes a 

periodic inflow of x with frequency / and initial (for t = 0) phase 0o. 0(t) 
is the Heaviside step function. We assume that the there is no inflow for 
t < 0, and it is switched on at t = 0. If 00 = 3 • 7r/2, then I[t) is a 
continuous function. In this case, /(t = 0) = 0. It then increases and finally 
oscillates. For any other phase, the inflow term is not continuous at t = 0; 
for example, if 0o = 7r/2 and then /(t = 0) = 2 ■ H. I[t) then decreases and 
hnally oscillates. The inflow term is always non-negative. For t > 0 the time 
average of I{t) equals A and is independent of the frequency and the initial 
phase. If the inflow amplitude A = 0, then (x = 0, ?/ = 0) is the only steady 
state of Fqs. fl2]l3l) and is stable. In the following analysis, we assume that the 
stable state of the system without flow is the initial state for the simulated 
evolution. 

Initially, let us consider the time evolution of the system for a constant 
inflow I{t) = A> 0 for t > 0 (thus, / = 0 and 0o = 0). The characteristics 
of the time evolution depend on the amplitude of the inflow term and on the 
initial state. In this case, the nullcline g{x,y,t) = 0 is the time-independent 
line with a definite slope determined by the value of a and a shift which 
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depends on the inflow amplitude A. Figured] shows the location of nullclines, 
calculated for e = 1, a = 0.55 and a few different values of the inflow. Let 
us assume that Ai and A 2 are the amplitudes for which the BN g{x,y) = 0 
nullcline is tangential to the sigmoidal-shaped nullcline h{x,y) = 0. The 
stable stationary states of the system can be located on two branches on 
the h{x,y) = 0 nullcline. One contains all of the points of the h{x,y) = 0 
nullcline located between point (0, 0) and the tangency point {xi,yi). We will 
call it the lower stable branch (LSB). The other is the upper stable branch 
(USB), and is formed by all of the points of the h{x,y) = 0 nullcline located 
above (x 2 , |/ 2 )- The stationary states located on the nullcline between {xi,yi) 
and (x 2 , 1 / 2 ) are unstable. In the case when A < A 2 , the only stationary state 
is located on the lower stable branch, so the system converges to the stable 
state yoo = hnit_j.oo 2/(i) such that yoo < yi regardless of the initial state. 
Similarly, for A > Ai, the single BN stationary state is located on the upper 
stable branch, and for all initial states the system converges to the stable 
state 1/00 > 1 / 2 • For A 2 < A < Ai, the stationary state that is approached 
for f > 00 depends on the initial state and on the partition of the phase 
space determined by the separatrices of the saddle point which is located 
on the middle branch of the h nullcline. This analysis also applies when 
the frequency of inflow oscillations is very high. In such a case, the flow 
oscillations are much faster than both the system dynamics and the system 
responses to the time-averaged value of the inflow A. 

For sufficiently slow oscillations of the inflow ( 0 < / <^ 1), the system 
can follow the slowly relocating stable state, the position of which varies 
according to the instantaneous value of the inflow. If the initial state of 
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the system is (a;(0) = 0,|/(0) = 0) and 2 ■ A < Ai, then y{t) < yi for all 
t. Therefore, the system state oscillates along the lower stable branch of 
the h{x, y) = 0 nullcline with the period dehned by the frequency of inflow 
oscillations. 2 ■ A > Ai, then there are intervals of time within which 
the system has a single stationary state located on the upper stable branch. 
During a single oscillation cycle, there are moments of time ti and ^2 at which 
y{ii) ^ Ui and 1 /(^ 2 ) > 2 / 2 , and thus oscillations that extend over both stable 
branches are expected. 

For moderate values of /, the system dynamics are too slow to closely 
follow the changes in the inflow value. In such a case, oscillations around a 
stable state located in the lower stable branch that extend to the unstable 
branch, as well as oscillations around a stable state located in the upper 
stable branch that extend to the unstable branch, should be observed. This 
is conhrmed by numerical simulations. 

The complexity of oscillations observed for a constant inflow amplitude 
A = 0.12 (thus A > Ai/2 but A < Ai ) as & function of flow frequency is 
illustrated in Fig. 2. As discussed above, for the selected amplitude and a 
low frequency of inflow oscillations, the system dynamics follow the time- 
dependent stationary state and oscillations of y{t) extend over both stable 
branches of the h{x, y) = Q nullcline. For intermediate frequencies, oscilla¬ 
tions accumulate on the upper stable branch of the nullcline and the min¬ 
imum value of y{t) increases with frequency. Then, at a certain frequency 
fc ( for the selected amplitude of oscillations, fc ~ 0.0312 ), the oscillations 
switch from the upper to the lower stable branch of the h(x, ?/) = 0 nullcline. 
The transition between oscillations located on different stable branches is 
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quite pronounced and should be easily detected in experiments with a sys¬ 
tem exhibiting hysteresis. Therefore, it becomes apparent that a cooperative 
system can discriminate the frequency of a perturbation if its amplitude 
remains hxed. Numerical simulations have also demonstrated that the fre¬ 
quency of the transition between oscillations on the USB and LSB depends 
on the phase 0o- The right upper corner of Fig. 2 shows two types of oscil¬ 
lations that are observed for / = 0.031. If 0o = 0, the system oscillates at 
the upper stable branch, but if 00 = ^, oscillations around the lower stable 
branch are seen. Fortunately for the application of this approach to discrim¬ 
ination, the interval of frequencies within which phase-dependent evolution 
is observed is very narrow. For A = 0.12, it is [0.0306, 0.0312]. The width of 
this interval (A/ ~ 0.0006) dehnes the precision in frequency discrimination. 

The dynamical system considered here can also be used to discriminate 
the amplitude of an applied perturbation. Figure 3 shows the time evolu¬ 
tion of y(t) for a few values of perturbation amplitude A and a hxed inhow 
frequency (/ = 0.05). As expected, for small amplitudes (A < 0.1122) the 
oscillations of y(t) are limited to the LSB. For larger amplitudes (0.1122 < 
A < 0.1361), the range of observed values of y(t) increases, but the oscilla¬ 
tions are still anchored on the LSB. Finally, if 0.13617 < A, the oscillations 
move onto the USB. The transition between the different types of oscillation 
is quite pronounced and can be used to discriminate the value of amplitude. 
Here, similar to the cases illustrated in Fig. 2, we observe a narrow interval 
of amplitudes ( AA ~ 0.0001) within which the type of oscillation depends 
on the initial phase. 


a = 0.55, 8 = 1 



Figure 1: Positions of nullclines for the dynamical system defined by 
Eqs. fl2]l3|) . The model parameters are: e = 1, a = 0.55, / = 0 and 
00 = 0. The nullcline h{x, ?/) = 0 is plotted with a solid line. The 
nullcline g{x,y,t) = 0 is shown for a few cases: A = 0 ( dotted line), 
A = A 2 = 0.02603 ( short-dashed line), A = Ai = 0.16445 ( long-dashed 
line). For the selected parameters of the model, the variables at tangential 
points are (xi,yi) = (0.47368,0.09607) and {x 2 ,y 2 ) = (1.23531,0.65339) 
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Figure 2: Time evolution of the dynamical system described by Eqs. fOIS]) as a 
function of the inflow frequency / for a hxed amplitude A = 0.12. The model 
parameters are: e = 1, a = 0.55. Tics and numbers on the frequency scale 
mark transitions between different types of oscillation. The initial phase 
is 00 = 37r/2 for all cases except the central one, for which 0o = 0. The 
horizontal green lines mark the values of yi and y 2 - 
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Figure 3: Time evolution of the dynamical system described by Eqs. fl^llS]) as 
a function of the inflow amplitude A for a hxed frequency / = 0.05. Tics 
and numbers on the frequency scale mark transitions between different types 
of oscillation. The initial phase is 0olT 37r/2 for all cases except at the top 
in the right column, for which (^q = 0. The horizontal dashed lines mark the 
values of yi and y 2 - 




























To give a more precise description of system evolution, let us introduce a 
classification of oscillations based on the minimum and maximum values of 
y{t) observed over a long time interval for which the evolution has reached a 
stationary state. We dehne: 


ymin — '^'^'^te[tmin,tmax] 


(4) 


and 


Umax V 


(5) 


Initially, we used tmin = 1000, where tmax = + 1000. Next we repeated 

the calculations for tmin = 2000. If there is a signihcant discrepancy in ymin, 
Umax obtained for these time intervals, then the procedure is repeated with 
tmin increased by an additional 1000 time units until agreement is attained. 

The type of oscillation is classihed through the comparison of ymin and 
Umax with the values of yi and 1 / 2 , as illustrated in Fig. 4. The classihcation 
of oscillations is summarized in Table I. 


Table I 


oscillation class 


condition 
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As shown in Fig. 4, the transitions between oscillation types (/) — 
\LSB -H- (/) —2 and (/) —2 -H- {III) are continnous because they result from 
an increase or decrease in the Umax value. Similarly, the transitions between 
oscillation types {HI) -H- {II) — 2 and {II) — 2 {II) — lUSB are conti- 
nous because they are related to an increase or decrease in the ymin value. 
In actual experiments, these transitions are difficult to detect because they 
require highly accurate data acquisition. On the other hand, the transition 
(/) — 2 -H- {II) — lUSB, on which the discrimination is based, can be easily 
detected because it is related to a discontinuous jump between Umax < y 2 
and ymin ^ y 2 - 
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Figure 4: Geometrical illustration of types of oscillation in the classification 
based on Umin and Umax- The dotted lines mark yi = 0.09607 and 1/2 = 
0.65338. 
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Figure 5 illustrates the regions of parameters (/, A) for which a given os¬ 
cillation pattern is observed in a model characterized by a = 0.55 and e = 1. 
The thick line separates the region of the phase space (/, A) in which class 
(I)-2 oscillations are observed, from the region where class (II)-l USB oscilla¬ 
tions appear. Let us denote points on this line as (/c, A^). The amplitude Ac, 
when treated as a function of /c, is a continuous, monotonically increasing 
function Ac = G{fc)- Therefore, the inverse function fc = G~^{Ac) exists. 
Our discrimination method is based on the determination of conditions in 
which a small change in fc or Ac qualitatively changes the character of the 
time evolution to force a transition between (J) — 2 and (//) — lUSB type 
oscillations. Let us assume that we want to measure the unknown inflow 
frequency and that we can regulate the inflow amplitude. The following 
procedure can be applied. Initially, we set a low amplitude so the system 
oscillates on the LSB (type {II) — 1 oscillations). Next, the amplitude is 
increased up to the moment A^ when oscillations of type {II) — lUSB are 
detected. The frequency of inflow can be estimated as fz = G~^{Az). This 
method works for all frequencies greater than /q, which corresponds to the 
tip of the {II) — lUSB region ((/o,Ao)). The accuracy of the estimation 
depends on the frequency and is high where the amplitude Ac is a rapidly 
increasing function of fc, here for 0.02 < fc S 0.1. This system can also be 
used to determine the amplitude of inflow when we can control the frequency. 
Now we set a low frequency and the system exhibits type {III) oscillations. 
Next, the frequency is increased up to the moment fy when oscillations of 
type (/) — 2 are detected. The amplitude of the inflow Ay is Ay = G{fy). 
Unlike for frequency, the range of discriminated amplitudes does not extend 
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outside the interval [Aq^Ai]. 

The phase diagrams, similar to that in Fig. 5 but for e = 1/5 and e = 5, 
are shown in Fig. 6 and Fig. 7 respectively. The results are qualitatively 
identical to those in Fig. 5, suggesting that the described changes in the 
system oscillations are generic and should also apply to other systems with 
hysteresis influenced by a periodic perturbation. 
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Figure 5: Phase diagram showing the type of oscillation as a function of 
inflow parameters {f,A). The horizontal lines indicate Ai = 0.16445 and 
Ai/2. The model parameters are e = 1 and a = 0.55 

. The thick solid line marks the boundary between oscillation classes (/) — 
2 and {II) — lUSB. The transition between these oscillations is used to 
determine the parameters of inflow. 
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Figure 6: Phase diagram showing the type of oscillation as a function of 
inflow parameters {f,A). The horizontal lines indicate Ai = 0.16445 and 
yli/2 . The model parameters are e = 1/5 and a = 0.55 
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Figure 7: Phase diagram showing the type of oscillation as a function of 
inflow parameters {f,A). The horizontal lines indicate Ai = 0.16445 and 
Axl‘1 . The model parameters are e = 5 and a = 0.55 
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3 Conclusions 


We have described how the time evolution of a cooperative system is de¬ 
pendent on the frequency and amplitude of a periodic stimulus. There is a 
narrow range of these parameters within which the characteristics of this evo¬ 
lution change in a qualitative manner: oscillations aronnd one stable branch 
change into oscillations on another branch. This phenomenon can be used to 
determine the amplitude or frequency of an applied perturbation. As for the 
numerical framework, we evaluated the effect of the rhythmicity of substrate 
input in a model biochemical system with sigmoidal kinetics, i.e. n = 3 in 
the Hill equation. Using numerical simulations, we separated the phase space 
of inflow parameters (amplitude and frequency) into regions where specihc 
types of oscillation are observed. The bonndary line separating oscillations 
with signihcantly different behaviors (type (/) — 2 and type (//) — lUSB 
oscillations) was identihed. The freqnency that causes a transition appears 
in a monotonic function of the inflow amplitude. The system can be used 
to determine the inflow frequency if we can control the inflow amplitude. It 
can also be used to determine the inflow amplitude when we can control the 
frequency. In other words, sigmoidal kinetics with the Hill equation can act 
as an inflow discriminator. 

This paper describes a system in which the nonlinear term in the kinetic 
equation for the y{t) variable is described by term (cf. Eq.(|3])) and the 
periodic inflow is described by a trigonometric function. We believe that 
these results are general, and qualitatively similar behavior can be expected 
in other systems with cooperative characteristics. We performed nnmerical 
simulations for a model based on Eqs.( [2]l3|l bnt with the inflow term in the 
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form J{t) = A ■ (tanh (7 • sin(27r/t + 0o)) + 1) ' 0(^) for different values of 
7 . Such periodic inflow becomes a square-like wave for large 7 . The phase 
diagrams that illustrate the type of oscillation as a function of / and A 
are qualitatively the same, as presented in Figs. (5-7). We also considered 
other nonlinear terms in the kinetic equation for y{t), like tanh (x — Xq) or 
1/(1 -I- exp (—5 • (x — xo)), and obtained similar results. Therefore, we be¬ 
lieve that real systems of chemical reactions with hysteresis can be used as 
discriminators in the manner described above. 

The present results can be regarded as a solution to the problem of the 
optimum stabilization of a system in an unstable state. Let us assume that 
Eqs. describe the time-dependent progress of a medical treatment where 

the variable y{t) represents the condition of a patient. The variable x{t) 
describes the time-dependent concentration of the curing drug. The states on 
the LSB and USB correspond to an ill and healthy patient, respectively. This 
simple model seems to realistically describe the basic features of drug therapy. 
It predicts that if the inflow of the drug is small, then the patient remains 
ill. Only a dose higher than a critical dose allows for successful treatment. 
However, some drugs are toxic ( such as those used in chemotherapy) and the 
total dose should be as small as possible. An analysis of the dynamical system 
presented in Fig. 5 can provide a solution: if we consider the periodic inflow 
of a drug in the form of Eq.(l), then the minimum amount of drug required 
to stabilize the patient in a healthy state corresponds to the bottom corner 
of the type (//) — lUSB oscillation region - here Aq = 0.1 and /o = 0.015. 
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